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Motivated by recent scanning tunneling microscopy experiments on surfaces of Bii-^Sb^/ [7, 10] 
and Bi2Te3,[8, 9] we theoretically study the electronic structure of a 3-dimensional (3D) topological 
insulator in the presence of a local impurity or a domain wall on its surface using a 3D lattice 
model. While the local density of states (LDOS) oscillates significantly in space at energies above 
the bulk gap, the oscillation due to the in-gap surface Dirac fermions are very weak. The extracted 
modulation wave number as a function of energy satisfies the Dirac dispersion for in-gap energies and 
follows the border of the bulk continuum above the bulk gap. We have also examined analytically 
the effects of the defects by using a pure Dirac fermion model for the surface states and found 
that the LDOS decays asymptotically faster at least by a factor of 1/r than that in normal metals, 
consistent with the results obtained from our lattice model. 



I. INTRODUCTION 

A three-dimensional (3D) topological insulator (TI) is 
a time-reversal invariant system with bulk gap but sup- 
port massless Dirac fermions with coupled spin and mo- 
mentum on the surface. [1, 2] The existence and the (odd) 
number of Dirac cones of the massless dispersion are pro- 
tected by the topology of the bulk band structure. [3] 
Recent first principle calculation reveals that Sb2Te3, 
Bi2Te3 and Bi2Se3 crystals are potential topological 
insulators. [4] Quite excitingly angle-resolved photoemis- 
sion (ARPES) experiments have confirmed many of the 
proposed topological insulators. [5] Before harvesting the 
novel properties of TPs, [6] an interesting question is how 
robust the surface Dirac fermions are against imperfec- 
tions, such as local and extended impurities. This issue 
is made interesting and challenging by recent scanning 
tunnelling microscopy (STM) measurements. [7-10] It is 
observed that while local impurities are seen to induce 
local density of states (LDOS) oscillation, [9] surprisingly 
the domain wall (a terrace on the surface) does not seem 
to induce LDOS oscillations for energies within the bulk 
gap. [7, 8, 10] 

We investigate the above issue by a model lattice 
hamiltonian for a 3-dimensional TI. The main results are 
as follows. 1) We obtain the surface Green's function 
and therefore the spectral function for the surface states. 
The contributions from surface Dirac fermions and the 
bulk extended states are clearly identified. 2) For an up- 
per terrace on the surface, we find the LDOS oscillation 
is vanishingly weak for energies below the bulk gap A, 
where Dirac fermions are best defined. The oscillation is 
significant above A, but it can be clearly ascribed to the 
contributions from the bulk extended states near the bot- 
tom of the conduction band. The numerical phenomenol- 
ogy is in nice agreement with the experiment. [7, 8, 10] On 
the other hand, the LDOS oscillation on the lower ter- 
race is globally weak. The asymmetry of the terraces can 
be attributed the the difference in the effective scatting 



mechanism. 3) For a local unitary impurity on a flat sur- 
face, we find LDOS oscillations at all energies, although 
it is relatively weaker below A. 4) Combining both ter- 
races and local impurities, we extract from the LDOS 
oscillations the energy uo dependence of the modulation 
wave number 2k u . For uo < A the dispersion (uo vs. kj) 
coincides with that of surface Dirac fermions, while for 
uo > A it is actually related to the bulk extended states. 
The above numerical results combine to support the ro- 
bustness of surface Dirac fermions. 5) For comparison, we 
substantiate analytical results using pure Dirac fermion 
models subject to impurity scattering. Asymptotically, 
the oscillation in LDOS, if present at all, has an energy 
dependent wave number 2k u , where k u is the on-shell 
momentum, and decays faster by a factor of 1/r than 
that for usual fermions. In particular, a hard domain 
wall in 2D Dirac models does not lead to any LDOS os- 
cillations at all. These are in qualitative agreement with 
the numerical results for the full 3d model. However, we 
do find and discuss differences between surface and pure 
Dirac fermions in connection to the nature of the wave 
functions. 

The rest of the paper is organized as follows. We dis- 
cuss the surface states using a 3D lattice model in sections 
II-IV, where a perfect surface, a terrace and a local impu- 
rity are discussed, respectively. Sec.V contains analytical 
results using pure Dirac models. We summarize and pro- 
vide remarks in connection to experiments in Sec. VI. 



II. SURFACE STATES IN A LATTICE MODEL 
OF TOPOLOGICAL INSULATORS 



We start with a 3D lattice model for topological insu- 
lators. The hamiltonian is given by H = ^kV^k^k, 
where ipk is a 4-spinor, k is the lattice momentum, and 

h k = [rn - ^ 2£ 6 (cos k b - 1)]T + ^ 2t b sin k b T b , (1) 
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FIG. 1: Schematic plot of a topological insulator with layered 
hexagonal structure, and a local impurity (filled circle) (a) 
and (upper and lower) terraces (b) on the surface. Here x, y 
and z denote the orthogonal axes, and the normal direction of 
the terrace edge is along x. The inter- layer hopping is along 
the vertical direction z. 

where m is a parameter controlling the bulk gap, % is 
the hopping amplitude along a bond b, fc& = k • b, and 
Tfr = (Tix + Y^y + r3z) • b. Here x, y and z are three 
orthogonal unit vectors, b is the unit vector along b, and 
T^'s (for \i = 0, 1, 2, 3) are Hermitian Dirac matrices that 
satisfy the Clifford algebra {r^,I\,} = 2S^ U . Explicitly 
we take Ti = <j\ <g) r 3 , T 2 = a 2 t 3 , T 3 = a r 2 and 
I\) = cr ri , where a and r are Pauli matrices (cr is the 
unit matrix). However, the general conclusion does not 
rely on the parameterization of the Dirac matrices. De- 
composing fok as /ik = ^2 £k,/iIV, one sees that the bulk 

dispersion is given by uo = ±^k with = \J^^ £k ^ 
We define A = mini^k as the bulk gap, which is assumed 
nonzero. We specify to a layered hexagonal lattice il- 
lustrated in Figs.l (the impurities are considered in the 
next section), as would be relevant to the experiments. [7- 
10] For simplest purposes we consider a nearest-neighbor 
tight-binding model and set 2t n = 1 and m = hence- 
forth. In this case A = 1. 

We use a recursion method to obtain the surface 
Green's function. In the presence of open surfaces normal 
to z, the planar momentum k|| = (k x , k y ) is still a good 
quantum number. For each k|| the hamiltonian can be 
decomposed formally into intra-layer part ^n^kj^n 

and inter-layer part ^ n V'n^V'n+i + h.c. where n is the 
layer index, 

fekfi = l m + 2t z - Yl 2 ^( cos h - i)]r + Yl 2tb sin kbTb > 

h^z h^z 

and h z = —t z (To + iTs) before setting our parametriza- 
tion. Define g = l/(zl — ) as the Green's function for 
an isolated layer (here z = uj + i0 + for retarded Green's 
function, and I is the 4x4 unit matrix), the surface 
Green's function for an N-layer lattice is given by, recur- 



sively, 

G^ 1 = g~ x - h\GN-\h z , (2) 

starting with G\ = g. Amazingly it is able to prepare 
samples with a series of layers in recent experiments. [11] 
In this case the recursion method is particularly useful to 
reveal the evolution of the surface states as the sample 
thickness increases. However, in this paper, we are only 
interested in infinite- layer lattices, for which G _1 = g~ x — 
h\Gh z holds for G = G^. 

The spectral function for the surface states is given by 

A(k||,^) = --ImTrG k|| , (3) 

where the planar momentum is indicated explicitly. The 
distribution in the momentum space is plot in Figs. 2(a)- 
(d) for a series of uj (only positive energies are displayed 
since the spectrum is particle- hole symmetric). It is seen 
that there is only a single dark spot at uo = 0. This 
indicates that we have but one Dirac cone in the Bril- 
louine zone. The spot evolves to a ring with increasing 
size for uo < A = 1. For uo > A bulk continuum starts to 
contribute, forming doubled ring structure at uo = 1.2 in 
combination with the Dirac ring. For still higher energies, 
the ring structure is blurred (not shown). Since A(k||,o;) 
is quite rotationally symmetric in the k|| space for the en- 
ergy window under concern (hexagonal structure shows 
up at higher energies), we plot it in Fig. 2(e) (gray) along 
a line cut (k x ,0). A clear massless Dirac dispersion is 
seen. The bulk continuum starts at uo = A and k|| = 0, 
and the outer border expands with increasing energies, 
while the Dirac dispersion persists for uo > A but even- 
tually diminishes for a; > 1.3. The contributions from 
surface Dirac fermions and bulk continuum are clearly 
separable. In comparison to pure Dirac models, the sur- 
face Dirac fermions have a momentum dependent spec- 
tral weight, the exact nature of which will be discussed 
elsewhere. 

III. TERRACES ON THE SURFACE 

Experimentally terraces often appear on the surface. 
We assume the upper terrace is one unit-cell higher than 
the lower one, as in the experiment [7-10] and illustrated 
in Fig. 1(b). The Green's function Q on the terraces is 
obtained as follows. First, we map the hexagonal lattice 
to a square lattice. The three principle translation axes 
in the hexagonal lattice are mapped to the horizontal, 
vertical and 45° axes on the square lattice, for which 
we denote the positions in real and momentum space 
as (u,v) and (p, g), respectively. Consider the top two 
layers. We denote 
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(a) oo=0 



(b) oo=0.4 




FIG. 2: (a)-(d): Intensity plot of A(k\\,cu) in the planar k||- 
space at oo = 0, 0.4, 0.8, 1.2. The intensity scales with the 
darkness. The view field is bounded by 2tt x 2tt (within the 
hexagonal Brillouine zone). Camera lighting is used to en- 
hance the contrast, (e) Intensity plot of A(k||,o;) along a cut 
(fcc,0) (gray). The open squares and circles are extracted 
from LDOS oscillations in Fig. 3(b) and Fig. 4, respectively. 
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FIG. 3: LDOS as a function of the normal displacement x on 
the lower (a) and upper (b) terrace, and as a function of the 
radial distance r from a local impurity (c). In each panel the 
energy increases from uj — at the bottom to uj — 2 at the 
top, with energy spacing Ao; = 0.05. The arrows indicate the 
lines at uj — 1 which coincides with the bulk gap. 



as the g-resolved unperturbed Green's function, where 
M is the number of lattice sites along the u- or v-axis, 
and a, (3 = U, L denotes the upper or lower layer. The 
Green's function G a @(p,q) can be obtained by coupling 
an isolated upper layer (through h z ) to a surface already 
considered in the previous section. Now setting m — > oo 
for u < on the highest layer effectively depletes the 
left half of that layer, leaving the lower-left (upper-right) 
layer mimic of the lower (upper) terrace. The perturbed 
Green's function Q for the top two layers can be obtained 
by the T- matrix formalism, [12] 

gf(u,u')=Gf(u-u') 

+ Gf(u-u a )T q (u a ,u b )G^(u b -u'), (5) 

U a ,b<0 

where the layer index U and the condition u a ^ < indi- 
cate that the the depleted region is the upper-left layer, 
and T q is the g-resolved T- matrix given by T~ 1 (u ai u b ) = 
-G^ u (u a - u h ). As what is required, Q^(u,u') = 
whenever ua or u'(3 falls on the upper- left layer. We em- 
phasize the Green's function obtained this way contains 
all effects from the bulk below the terraces. The LDOS 
on the terraces is given by, 

p a (u,u) = -^lmTr^g™(u,u), (6) 
q 

where a = L (a = U) for u < (u > 0). Finally we map 
u back to the normal displacement x from the edge on 
the hexagonal lattice, as shown in Fig. 1(b). 

The LDOS as a function of x is shown in Fig. 3(a) and 
(b) for the lower and upper terraces, respectively. In each 



panel, the energy increases from co = at the bottom up 
to u = 2 on the top, with energy interval Auj = 0.05. The 
arrow highlights the line with uo = A. For reasons to be 
clarified, we denned the origin of x differently for the two 
terraces. We observe that the oscillation of LDOS on 
the lower terrace is much weaker than that on the upper 
terrace. This important asymmetry can be checked ex- 
perimentally, and can be understood as follows. As far as 
the LDOS (or the local Green's function) is concerned, 
the upper terrace can be obtained by setting m^ooon 
the left half of a surface. Therefore the upper terrace ef- 
fectively experience a hard wall on the left. On the other 
hand, the lower terrace can be obtained by coupling a 
surface to an isolated upper-right layer through h z . Thus 
the lower terrace effectively experiences a soft boundary 
on the right. The asymmetry clearly follows from the dif- 
ference in the effective scattering mechanism, and makes 
it more sensible to define the origin of x as the boundary 
of such scattering in Figs. 3 (a) and (b). 

Let us concentrate on Fig. 3(b). Here the LDOS os- 
cillation is negligibly small for uo < A, while more and 
more peaks appear significant for uo > A. The LDOS 
oscillation is a manifestation of quasiparticle scatter- 
ing interference, [13] loosely referred to as the Friedel 
oscillation. [14] The wavelength A^, as is easily extracted 
from the first peak, should correspond to the wave num- 
ber 2k u , the characteristic momentum transfer during 
elastic quasiparticle scattering. Therefore = tt/X^. 
By this means we obtain a dispersion uo vs. fc^, which 
we plot in Fig. 2(e) (open squares). It clearly traces the 
outer border of the bulk continuum, and is therefore not 
related to the surface Dirac fermions. The lack of visi- 
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ble oscillations for uj < A is in nice agreement with the 
experiment. [7, 8, 10] The oscillation turns out to be more 
visible around a local impurity which we discuss in the 
next section. 



co=0.2 





co=0.4 




IV. A LOCAL IMPURITY ON THE SURFACE 

We consider a scalar impurity at the origin (r = 0) on 
the surface with the potential matrix VI, as illustrated 
in Fig. 1(a). (Non-scalar impurities [15] can be discussed 
along the same line but we assume scalar ones are more 
likely to occur.) Given the unperturbed surface Green's 
function G, the perturbed Green's function Q (in real 
space on the surface) can be obtained again by the T- 
matrix formalism, 



oo=0.6 



oo=0.8 



oo=1.0 



£(r, r') = G(r - r') + G(r)TG(-r'), 



where 



G(r) 



1 



5> 



C||T 



(7) 



(8) 



and T is given by T 1 
given by 

p{uj,v) = - 



V 1 ! - G(0). The LDOS is 



--ImTr£(r,r). 

7T 



(9) 



We set V — > oo for a unitary impurity. The LDOS along 
a principle translation axis is plot in Fig. 3(c). We see 
LDOS oscillation for uj > A similar to that in (b), al- 
though it is slightly weaker and has more complicated 
patterns. The oscillation below A is still weak but the 
peaks are visible for uj > 0.5 and they shift toward the 
origin with increasing energy. To have a better idea of 
the wavelength, we plot the LDOS map in Figs. 4 for a 
few values of uj. The view field bounded by 26 x 26 is 
extracted from a 400 x 400 surface. The oscillation, even 
if it is present, is beyond the view field for uj < 0.4, while 
it is clear (from the dark rings) for higher energies. We 
extract the wave number 2k UJ from the radii of the rings 
and plot uj vs. k^ in Fig. 2(e) (open circles) for uj < A. It 
clearly follows the Dirac dispersion. For uj > A, the os- 
cillation pattern is more complicated. In particular there 
are supermodulations for uj = 1.2 and uj = 1.4, which can 
be understood from the double-ring structure in the cor- 
responding map plot of the spectral function in Fig. 2(d). 



V. ANALYTICAL RESULTS FOR PURE DIRAC 
MODELS 

ID case: For illustrative purpose and to set up nota- 
tions, we start with the single-particle hamiltonian in a 
ID Dirac model, in the momentum space, h = /c<r n , where 






FIG. 4: Two-dimensional maps of the DOS near a local 
unitary impurity for specific values of uj. The view field is 
bounded by 26 x 26. The strength scales with the darkness. 
The DOS at the origin vanishes in our case, but it is arti- 
ficially reset to the value far from the impurity in order to 
enhance the contrast. Camera lighting is used. 



<j n is one of the three Pauli matrices and the Dirac ve- 
locity is set to unity henceforth. The unperturbed Mat- 
subara Green's function in the real space is given by 



G(x) 



I 



dk iuj n (jQ + k(j n 



ikx 



2tt loI + k 2 
-i [u! n a + |w n |sign(a;)] <f>(x), 



(10) 



where do is the unitary matrix, uj n is the Matsubara fre- 
quency, and we defined a kernel function 



(j)(x) 



I 



dk 



n ikx 



1 



2-xujl + k 2 2\uj n 



-\uj n x\ 



(11) 



The on-site Green's function G(0) = — iuj n (j)(0)ao turns 
out to be a scalar. Suppose there is a local impurity po- 
tential V(Jo at the origin, the perturbed Green's function 
is conveniently obtained by the T-matrix formalism, 



G(0) + G(x)TG(-x / ), 



(12) 



G{x,x') = 

where T _1 = V~ 1 cfq — G(0) ex do is the inverse of the 
T-matrix. We will concentrate on 5g(x), the trace of 
the change of the on-site Green's function as a function 
of x. Upon analytical continuation iuj n — > uj + i0 + , its 
imaginary part gives the change of the LDOS. For i/O, 
we find 



5g{x) = Tt[G(x)TG(-x)] 

oc Ti[(iuj n cro + i\uj n \a n )(iuj n a - i\uj n \a n )] 



= 0. 



(13) 



Therefore right away from the impurity site, the LDOS 
is unaffected at all. Moreover, using the above T- matrix 
formalism it is easy to verify that the off-site Green's 
function Q(x,x') = —G(x — x') for xx' < in the limit 
V — > oo. This signifies perfect transmission (albeit with 
a phase lag of 7r), a manifestation of the Klein Paradox. 
The mechanism behind this effect is the chirality (the 
alignment between the momentum and spin polarization) 
of the unperturbed eigen states. If an incoming state 
were scattered backward, energy conservation requires a 
flipping of the spin. But a scalar impurity can not flip 
the spin, so the scattering matrix element between these 
states vanishes identically. 

For comparison, the retarded Green's function in ID 
metal reads G(x) ~ —ie^^/v^ (k^ and v u are the on- 
shell momentum and group velocity at real energy uj). 
The oscillation can be worked out by the T-matrix ap- 
proach again. It has a wave number 2k u and does not 
decay at all in the clean limit implicitly assumed. (In 
reality, the oscillation must decay beyond the mean free 
path.) 

2D case: We write h = k • a = pcr x + qa y for 2D Dirac 
fermions. We first consider a domain wall with scalar po- 
tential V<jo along the y direction. Then q is still a good 
quantum number. The scattering geometry is illustrated 
in Fig. 5, where the circle of radius k u is an energy shell, 
radial arrows indicate the momentum dependent spin po- 
larization and the long thick arrow indicate a scattering 
process with incident angle (so that q — k u sin#). The 
(/-resolved unperturbed Green's function is given by 



G q (x) = - j 



dk iuo n (jQ + q<Jy + pa n 



2tt co 2 + p 2 + q 2 
ico n a + qcr y + iy/q 2 + u%(T x sign(x) <j> q {x\\4) 



where we defined 

<j> q {x) = 



-y/?+U*\x\ 



In particular, 

G q (0) = -(iu! n a + q<Jy)<f>q(0), 



(15) 



(16) 



is no longer a scalar unless q = 0. The (/-resolved per- 
turbed Green's function is given by [12] 



G q (x, x f ) = G q (x - x f ) + G q {x)T q G q {-x f ) 

with 

T" 1 = y-V - G ? (0). 

By straightforward algebra, we find 

8g q {x) = Tx[G q (x)T q G q (-x)} 
Aq 2 V~ l 



(17) 
(18) 



(19) 




FIG. 5: Schematic plot of an energy shell at uj with an on- shell 
momentum k^. The radial arrows indicate the spin polariza- 
tions. The long thick arrow indicates a scattering process 
that conserves q in the case of a domainwall with its normal 
direction along the p-axis, and is the incident angle. 



where we defined V q = V0 g (O) = V/2^q 2 + 

Clearly the factor of q 2 in Sg q (x) follows from the re- 
quirement of spin overlap sin 2 during elastic scattering, 
as illustrated in Fig. (5). This leads to Sg q (x) = for 
q = (normal incidence), recovering the ID Klein para- 
dox, as would have been anticipated. More importantly, 
Sg q (x) — > in the unitary limit V — > oo (a hard wall) 
for all values of q, and consequently in the real space 
Sg(x) = J Sg q (x)dq/27r — > as long as \x\ > 0. This ex- 
plains nicely the numerical results in Sec. Ill for the ter- 
race that the LDOS barely oscillates below the bulk gap 
energy. (In Sec. Ill we considered m — > oo, the effect of 
which is however identical to V — > oo here.) In the mean 
time, the off-site Green's function G q (x, x') = —G q (x—x') 
for xx' < holds, signaling perfect transmission in ana- 
logue to the ID case. This surprising result indicates that 
the Klein paradox also works in 2D but only for a hard 
wall. 

For general strength of V (or soft walls), the asymp- 
totic behavior in the limit of \uj n x\ ^> 1 (where a saddle 
point approximation is valid) is given by 



Sg(x) 



I 



Sg q (x) 
V 



(2 + iVsign(j n ) 2 7r\x\ y \^ n x 



-2\cj n x\ 



,(20) 



where q = k u sin is used. Upon analytical continuation 
\uo n \ — > iuo, this amounts to an oscillation of LDOS with 
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an energy-dependent wave number 2k^ = 2\co\ and an 
envelope function that decays as |x| -3 / 2 . 

For comparison, in a 2D metal, G q (x) ~ 
—ie lk "\ x \ cosG /v (Jl jCOs6 in real frequency. By the T- 
matrix formalism the domain wall leads to Sg(x) oc 
- $ dQe 2ik «\ x \ cosd /(V^VuCosQ + i) oc e 2ik -\ x \ / ^KM 
for k u \x\ ^> 1. So the LDOS oscillation decays as 
|&u;^| -1 / 2 near a domain wall in a 2D metal. [16] Notice 
that kuj = o is just the Fermi wave number. 

For a point impurity in 2D, the unperturbed Green's 
function for Dirac fermions is given by 



G(r) =- 



d 2 k iuj n cro + k • a 



ikr 



where we find 



0(r) 



(2tt) 2 u\ 
co n a + (\co n 

-\u n \r 



k 2 



— )r • a 

2r J 



47T 



2tt 



\Vn\r 



0(r), (21) 



(22) 



in the asymptotic limit \co n \r ^> 1. For the on-site 
Green's function, a cutoff A in momentum has to be in- 
troduced so that 



= -^o ln A 

47T CO, 



2 

2 ' 



(23) 



Using the T-matrix formalism with T 1 = V 1 ao — G(0), 
we find 



TrG(r)TG(-r) 
1 



-2|w„|r 



47TF" 1 +iw n ln(A 2 /w2) r 2 



(24) 



to leading order in l/\co n \r. Upon analytical continua- 
tion, we see that the change of LDOS oscillates with an 
energy dependent wave number 2A^ = 2\uj\ and decays as 
r -2 . Notice that the oscillation exists even for V — > oo, 
in contrast to the case of the domain wall. This is con- 
sistent with the results in Sec. IV. However, there is an 
important difference. Here a sharp resonance state ap- 
pears at co = (in the unitary limit V — > oo) as seen 
from analytical continuation of 5g(r), a well known re- 
sult in many contexts. [17] This resonance state is absent 
in Sec. IV, as can be seen in Fig. 3(b) or Figs. 4. Qualita- 
tively this is due to the fact that surface Dirac fermions 
are not completely confined on the surface. Depending 
on the planar momentum the wave function has varying 
extent of amplitude on the surface. Therefore the impu- 
rity is only partially seen by the surface Dirac fermions. 
Moreover, they exist only in a limited regime in the mo- 
mentum space. An accurate discussion of these issues are 
left for future studies. 

For comparison, in a 2D metal, G(r) = 



-ifdOe 



iku,r cos 



k w /27TV u oc e <fe - r (fe w r)- 1 / 2 for k u r > 1. 



By the T-matrix approach a local impurity leads to 
6g(r) oc G(r) 2 oc e 2ik « r /k^r for k u r > 1. So the LDOS 
oscillation decays as r _1 near a local impurity in a 2D 
metal. [16] 

Similar analysis could be proceeded for 3D Dirac 
fermions in the presence of 2-, 1- and 0-dimensional impu- 
rities. As the impurity states would not be easily probed 
by STM, we do not go into detailed analysis. Some in- 
teresting behaviors are as follows. First, a 2D hard wall 
does not lead to any oscillation in LDOS away from the 
wall. For lower-dimension or general scalar impurities 
the oscillation is present, and decays faster than that for 
normal metals by a factor of 1/r. 



VI. SUMMARY AND REMARKS 

We analyze the behavior of surface states of topological 
insulators against local and extended impurities. Using 
an effective 3D lattice model we first obtain the spectral 
function of the surface states. This enables us to identify 
the contributions from massless Dirac fermions and bulk 
continuum. We then study the effects of a local impurity 
and a terrace on the surface. Away from a local impurity, 
the spatial oscillation in LDOS below the bulk gap A is 
much weaker than that above A. The LDOS oscillation 
is globally weak on the lower terrace. While on the upper 
terrace, the LDOS oscillation is barely visible below A, in 
perfect agreement with the STM measurement, but it is 
significant above A. The asymmetry of the LDOS oscilla- 
tion on the lower and upper terraces is attributed to the 
difference in the effective scattering mechanism. From 
the LDOS oscillations we extract the modulation wave 
number 2k UJ as a function of energy co. The dispersion 
(co vs. k u ) follows that of Dirac fermions for co < A, but 
it follows the border of the bulk continuum for co > A. 
The numerical results combine to reveal that the surface 
Dirac fermions are rather immune to the imperfections. 
We discuss such a behavior analytically by pure Dirac 
models. Because of a cancellation due to the alignment 
of momentum and spin polarization, the 2k u oscillation 
in LDOS, if present, decays asymptotically faster by a 
factor of 1/r than that for usual fermions. Specifically, 
the oscillation is absent for a hard domain wall. Such 
behaviors are consistent with the 3d lattice model. We 
also find and discuss differences between surface and pure 
Dirac fermions. 

Before closing, some remarks are in order. First, 
our simple model does not yield wavy energy shells (in 
the momentum space) for surface Dirac fermions due to 
warping terms, which would arise if longer-range hop- 
ping are considered. The warping term produces par- 
tial nesting and is responsible for fine features in Bi2Te3. 
However, the merit of our simple model is it uncovers 
material-independent qualitative difference between sur- 
face Dirac fermions and bulk continuum. A material- 
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dependent calculation will be provided elsewhere. Sec- 
ond, we notice that a real terrace edge is not likely a 
straight line. This may be mapped to a straight edge but 
with excess impurities nearby. According to our results of 
both local impurity and terrace, a wiggling terrace edge 
should also lead to visible oscillations below the bulk gap, 
consistent with experiments. [9] Finally, in a model with 
multiple Dirac cones for the surface states, elastic scatter- 
ing between cones is not suppressed by the chirality, and 
is therefore expected to induce visible LDOS oscillations 
even below the bulk gap. 

While finalizing the writing of this work, we become 
aware of a related work using a 2D continuum model 
with warping terms. [18] 
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